DALS020-距离与降维

前言

这一部分是《Data Analysis for the life sciences》的第8章统计模型的第1小节，这一部分的主要内容涉及降维分析的一些原理，例如SVD，投影，旋转等，相应的Rmarkdown文档可以参考作者的Github。

距离(distance)的概念非常直接，例如，当我们把动物聚为亚群时，我们其实就是隐含地定义了一个距离，从而使我们可以说亚群内的动物彼此“接近”，如下所示：

我们使用的分析高通量数据的方法地都直接或间接地与距离有关。许多聚类方法和机器学习方法都需要使用特征值或预测因子来定义距离。例如热图是基因组学与高通量数据领域里使用最为广泛的工具，如果我们要生成热图，就需要明确计算距离，如下所示：

在上面的热图中，每个方格代表的数值储存在一个矩阵里，每具方格的行与列被聚类后（注：热图常用红色与绿色表示，但是对色盲人士来说，这两种颜色是最难分辨的颜色）用不同的颜色表示。在这一部分中，我们将学习必要的数学知识与计算技能来了解和创建热图。我们先来回顾一下数学上对距离的定义。

欧氏距离(Euclidean Distance)

现在我们在一个笛卡尔坐标系(Cartesian plane)中定义A点与B点的距离，如下所示：

library(rafalib)
mypar()
plot(c(0,1,1),c(0,0,1),pch=16,cex=2,xaxt="n",yaxt="n",xlab="",ylab="",bty="n",xlim=c(-0.25,1.25),ylim=c(-0.25,1.25))
lines(c(0,1,1,0),c(0,0,1,0))
text(0,.2,expression(paste('(A'[x]*',A'[y]*')')),cex=1.5)
text(1,1.2,expression(paste('(B'[x]*',B'[y]*')')),cex=1.5)
text(-0.1,0,"A",cex=2)
text(1.1,1,"B",cex=2)

其中，欧氏距离的定义如下所示：

$$\sqrt{ (A_x-B_x)^2 + (A_y-B_y)^2}$$

高维数据的距离

现在我们使用一个数据集，这个数据集中含有189个样本，22215个基因，如下所示：

library(devtools)
# install_github("genomicsclass/tissuesGeneExpression")

library(tissuesGeneExpression)
data(tissuesGeneExpression)
dim(e) ##e contains the expression data
table(tissue) ##tissue[i] tells us what tissue is represented by e[,i]

这些数据代表了8个组织（每个组织中有多个样本）的RNA表达水平，如下所示：

> dim(e) ##e contains the expression data
[1] 22215   189
> table(tissue) ##tissue[i] tells us what tissue is represented by e[,i]
tissue
 cerebellum       colon endometrium hippocampus      kidney       liver    placenta 
         38          34          15          31          39          26           6

我们现在描述一下在这个数据集中，不同样本之间的距离。我们也许对在不同样本中表达相似的基因感兴趣。

为了定义这个距离，我们需要知道这些点是什么，因为我们计算数学上的距离需要这些点。由于这个数据集是高维数据集，这些点就无法直接放在笛卡尔坐标系中。相反，我们会把这个数据集放在更高维度的坐标系中。例如样本$i$是由22215维空间的一个点定义的（这个空间可以写为$(Y_{1,i},\dots,Y_{22215,i})^\top$）。特征值$g$是由一个189维空间的一个点定义的（这个空间可以写为$(Y_{g,1},\dots,Y_{g,189})^\top$）。

一旦我们定义好了这些点，那么欧氏距离就可以使用我们前面类似的方法进行计算，例如，两个样本$i$和$j$的距离为：

$$\mbox{dist}(i,j) = \sqrt{ \sum_{g=1}^{22215} (Y_{g,i}-Y_{g,j })^2 }$$

两个特征值$h$和$g$的距离为：

$$\mbox{dist}(h,g) = \sqrt{ \sum_{i=1}^{189} (Y_{h,i}-Y_{g,i})^2 }$$

矩阵代数与距离

样本$i$和$j$之间的距离可以写为：

$$\mbox{dist}(i,j) = (\mathbf{Y}_i - \mathbf{Y}_j)^\top(\mathbf{Y}_i - \mathbf{Y}_j)$$

其中，$\mbox{Y}_{i}$和$\mbox{Y}_{j}$代表第$i$列和第$j$列。这种写法在实际计算中非常方便。

案例

现在我们使用上面的矩阵代数来计算一下距离。现在我们计算样本1与样本2（它们都是肾脏组织）的距离，然后再过计算样本87的距离（结肠），如下所示：

x <- e[,1]
y <- e[,2]
z <- e[,87]
sqrt(sum((x-y)^2)) # Kindey
sqrt(sum((x-z)^2)) # Colon

结果如下所示：

> sqrt(sum((x-y)^2)) # Kindey
[1] 85.8546
> sqrt(sum((x-z)^2)) # Colon
[1] 122.8919

从结果中我们可以发现，肾脏组织之间距离比较近。另外一种计算距离更快的方式就是使用矩阵代数，如下所示：

sqrt( crossprod(x-y) )
sqrt( crossprod(x-z) )

结果如下所示：

> sqrt( crossprod(x-y) )
        [,1]
[1,] 85.8546
> sqrt( crossprod(x-z) )
         [,1]
[1,] 122.8919

现在我们一次计算出所有的距离，我们可以使用dist()函数，这个函数可以计算出每行之间的距离，现在我们感兴趣的是不同样本之间的相似性，因此我们需要使用t()来转换一下矩阵，如下所示：

d <- dist(t(e))
class(d)

结果如下所示：

> d <- dist(t(e))
> class(d)
[1] "dist"

从结果中我们可以发现，计算结果d是一个dist类，为了得到它的具体数值，我们需要将其强行转换为矩阵，并使用索引操作，如下所示：

as.matrix(d)[1,2]
as.matrix(d)[1,87]

结果如下所示：

> as.matrix(d)[1,2]
[1] 85.8546
> as.matrix(d)[1,87]
[1] 122.8919

这里我们需要注意的是，我们对数据集e使用了函数dist()，这个函数计算的是基因之间所有两两矩离，它最终会形成一个$22215\times 22215$的矩阵。

练习

P322

降维操作

相应的Rmarkdown见作者的Github。

可视化数据是分析高通量数据中最重要的步骤之一。正确的可视化方法可能会发现实验数据的问题，这些数据可以呈现标准分析的结果。我们已经展示了可视化数据的全局方法，但是由于数据的高维特性，使得发现列之间或行之间关系的图形变得复杂。例如，如果要比较189个样本之间的特性，我们必不得不创建17766个MA图。创建一个单独的散点图明显不合适，因为数据量太大。

我们将介绍其于降维的探索性数据分析的强大技巧。一般的想法就是将数据集降至较低维度的同时又保留重要的特性，例如样本之间的距离。如果我们能够将数据降低到2维，那么我们就能很容易地画出图形。降维的背后就是奇异值分解(Singular value decomposition, SVD)，这种思路也可以应用于其他情况。在介绍SVD背后的复杂数学原理之前，我们将会使用一个简单的案例来介绍一下它的思路。

案例：将2维数据降低至1维

现在我们来看一个案例，这个案例是有关双胞胎身高的。我们来模拟生成100个二维数据点，它们表示每个人与其均值的偏离的标准差的数目，每对数据点表示一对又胞胎，如下所示：

library(rafalib)
library(MASS)
set.seed(1)
n <- 100
y=t(mvrnorm(n,c(0,0), matrix(c(1,0.95,0.95,1),2,2)))
mypar()
plot(y[1,], y[2,], xlab="Twin 1 (standardized height)", 
     ylab="Twin 2 (standardized height)", xlim=c(-3,3), ylim=c(-3,3))
points(y[1,1:2], y[2,1:2], col=2, pch=16)

为了辅助说明问题，我们可以将上面的模拟数据视为高通量基因表达数据，其中双胞胎的配对数据表示了N个样本，双胞胎的2个身高表示基因表达数据。我们是对任意2个之间的距离感兴趣。我们可以使用dist()函数来进行计算。例如，上图的2个橙色数据点的距离为：

d=dist(t(y))
as.matrix(d)[1,2]

结果如下所示：

> d=dist(t(y))
> as.matrix(d)[1,2]
[1] 1.140897

如果两维数据太复杂（这里只是假设），我们只想制备一维图，那怎么办呢？例如，我们能否将这些数据减化为一维矩阵，同时保留这些点与点之间的距离信息呢？

如果我们回顾再来看这张图，在任何一对数据点之间画一条线，那么这条线的长度就是这两点这之间的距离。这些线倾向于沿着对角线的方向分布。我们以前到过MA图，这种图就是将原始散点图的对角线“旋转”了一下，将原来的对角线旋转到与x轴平行的位置形成的，如下所示：

z1 = (y[1,]+y[2,])/2 #the sum
z2 = (y[1,]-y[2,]) #the difference
z = rbind( z1, z2) #matrix now same dimensions as y
thelim <- c(-3,3)
mypar(1,2)
plot(y[1,],y[2,],xlab="Twin 1 (standardized height)",ylab="Twin 2 (standardized \
height)",xlim=thelim,ylim=thelim)
points(y[1,1:2],y[2,1:2],col=2,pch=16)
plot(z[1,],z[2,],xlim=thelim,ylim=thelim,xlab="Average height",ylab="Differnece \
in height")
points(z[1,1:2],z[2,1:2],col=2,pch=16)

上图的左图就是原始的散点图，右图则是MA图。

在后面章节中，我们将会使用线性代数来表示这些数据的变换（也就是旋转）。这里我们可以通过将y的相乘来得到z，如下所示：

$$A = \, \begin{pmatrix} 1/2&1/2\\ 1&-1\\ \end{pmatrix} \implies z = A y$$

我们将两侧都乘以$A^{-1}$，则得到z，如下所示：

$$A^{-1} = \, \begin{pmatrix} 1&1/2\\ 1&-1/2\\ \end{pmatrix} \implies y = A^{-1} z$$

旋转

在上图相，相对于其它点之间的距离，两个橙色上炽之间的距离大致保持一致。所以的点其实都是如此。对上面转换进行简单的重新缩放，将会使前后的距离完全相同。我们要做的就是乘以一个标量，从而保留每个数据点的标准差。如果你认为y的列是一个独立随机变量，其标准差为$\sigma$，那么我们要注意到$M$与$A$的标准差如下所示：

$$\mbox{sd}[ Z_1 ] = \mbox{sd}[ (Y_1 + Y_2) / 2 ] = \frac{1}{\sqrt{2}} \sigma \mbox{ and } \mbox{sd}[ Z_2] = \mbox{sd}[ Y_1 - Y_2 ] = {\sqrt{2}} \sigma$$

这就说明，如果我们将上面的转换变为如下形式：

$$A = \frac{1}{\sqrt{2}} \begin{pmatrix} 1&1\\ 1&-1\\ \end{pmatrix}$$

那么$Y$列的SD就会变得与$Z$列的方差一样。此外，我们要注意到，$A^{-1}A=I$。我们称这种特性为正交(orthogonal)，并且它保留了上述SD的特性。因此就保留了距离信息：

A <- 1/sqrt(2)*matrix(c(1,1,1,-1),2,2)
z <- A%*%y
d <- dist(t(y))
d2 <- dist(t(z))
mypar(1,1)
plot(as.numeric(d),as.numeric(d2)) #as.numeric turns distnaces into long vector
abline(0,1,col=2)

我们称这种转换为y的旋转：

mypar(1,2)
thelim <- c(-3,3)
plot(y[1,],y[2,],xlab="Twin 1 (standardized height)",ylab="Twin 2 (standardized height)",xlim=thelim,ylim=thelim)
points(y[1,1:2],y[2,1:2],col=2,pch=16)
plot(z[1,],z[2,],xlim=thelim,ylim=thelim,xlab="Average height",ylab="Differnece in height")
points(z[1,1:2],z[2,1:2],col=2,pch=16)

我们之所以优先使用这种转换，是因为我们注意到所有的点是沿着对角线进行分布的，我们将对角线进行转换后，对角线与x轴平行。所以这个旋转实际上就达到了我们最初的要求：我们只需要一个维度就可以保留点与点之间的距离。现在让我们删除了第二个维度z，并重新计算一下距离：

d3 = dist(z[1,]) ##distance computed using just first dimension
mypar(1,1)
plot(as.numeric(d),as.numeric(d3))
abline(0,1)

仅用一维数据进行的距离计算就很接近实际距离，并且降低了维度，将2维降低到了1维。转换后的数据的第1约就是第一主成分。这一思想促进了使用主成分分析(PCA)和奇异值分解(SVD)来实现更广泛的降维。

关于与其他解释区别的重要说明

如果你在网上搜索PCA的描述，你会注意到网上的描述与我们这里描述在符号上有些出入。这主要是因为在PCA中，通常使用行来表示实验单元（就是样本）。因此，在我们这里的实验中，$Y$通常会被转换为$N \times 2$矩阵。在统计学中，这也是最为普遍的表示数据的方式：每行表示一个样本。然而，由于实际原因，在遗传学中，通常使用列表示样本。例如行表示基因，列表示样本。由于这个原因，在这本书中，我们会解释PCA和与之相应的数学计算会与常规的方式有所不同。因此，在网上找到的相关的许多PCA的解释都是先从样本的协方差矩阵开始的，它通常使用$\mathbf{X}^\top\mathbf{X}$表示，并且每个单元格表示两个实验单元之间的协方差。然而，要做到这一点，我们需要使用$\mbox{X}$的行表示实验单元。因此，在我们上面的符号中，在经过缩放后，你必须要使用$\mathbf{Y}\mathbf{Y}^\top$来进行计算。总之，如果你想让我们的解释与其他有关的PCA内容相符，就必须对这本书中使用的矩阵进行转置。

奇异值分解

相关的Rmarkdown参考作者的Github。

在前面的部分中，我们展示了降维分析，以及如何让我们使用一维数据来替代原来的二维数据，近似地表示点与点之间的距离。奇异值分解(SVD)是我们这种方法的推广。在这个案例中，SVD对原始数据进行了转换。这种转换具有一些非常有用的属性。

SVD计算的主要结果就是，我们可以写为一个$m\times n$矩阵，对于矩阵$\mbox{Y}$就写为：

$$\mathbf{U}^\top\mathbf{Y} = \mathbf{DV}^\top$$

其中，

• $\mathbf{U}$ 是一个 $m \times p$ 正交矩阵；
• $\mathbf{V}$ 是一个 $p \times p$ 正交矩阵；
• $\mathbf{D}$ 是一个 $n \times p$ 对角矩阵。

其中，$p=\mbox{min}(m,n)$，$\mathbf{U}^\top$对数据$\mathbf{Y}$进行旋转，这一步非常有用，因为$\mathbf{U}^\top \mathbf{Y}=\mathbf{VD}$列的变异（精确的平方和）会下降。因此$\mathbf{U}$是正交辞职，我们可以将SVD写为如下形式：

$$\mathbf{Y} = \mathbf{UDV}^\top$$

实际上这个公式更加普遍，我们也可以将转换写为如下形式：

$$\mathbf{YV} = \mathbf{UD}$$

$Y$的这种转换也会生成一个矩阵，这个矩阵的列的乘方和是递减的。

将SVD应用到我们的案例中，就是如下结果：

library(rafalib)
library(MASS)
n <- 100
y <- t(mvrnorm(n,c(0,0), matrix(c(1,0.95,0.95,1),2,2)))
s <- svd(y)

我们可以马上就是看到使用了SVD后生成的转换后的矩阵非常类似于我们前面的案例中的结果：

round(sqrt(2) * s$u , 3)

结果如下所示：

> round(sqrt(2) * s$u , 3)
       [,1]   [,2]
[1,] -0.982 -1.017
[2,] -1.017  0.982

当我们旋转后，绘制成的图形称为主成分(principal coimponent)：这里只绘制出了第一个主成分和第二个主成分。如果我们想要从SVD中获取主成分，只需要旋转后的$\mathbf{U}^\top\mathbf{Y}$ 即可：

PC1 = s$d[1]*s$v[,1]
PC2 = s$d[2]*s$v[,2]
plot(PC1,PC2,xlim=c(-3,3),ylim=c(-3,3))

用途

使用SVD的用途并不十分明显，我们可以看一些案例。在这个案例中，我们将会极大地降低$V$的组倒数，并且仍然能够构建$Y$。

现在我们来对基因表达谱进行SVD的计算，我们可以只使用表达谱中的100个基因的子集，这样计算会快一点，如下所示：

library(tissuesGeneExpression)
data(tissuesGeneExpression)
set.seed(1)
ind <- sample(nrow(e),500)
Y <- t(apply(e[ind,],1,scale)) #standardize data for illustration

使用svd()函数可以返回3个矩阵（D矩阵仅返回对角线元素），如下所示：

s <- svd(Y)
U <- s$u
V <- s$v
D <- diag(s$d) ##turn it into a matrix

我们首选要注意到，我们可以重构$y$：

Yhat <- U %*% D %*% t(V)
resid <- Y - Yhat
max(abs(resid))

结果如下所示：

> max(abs(resid))
[1] 3.552714e-14

如果我们看一下$\mathbf{UD}$的平方和，我们会看到最后几个非常接近于0（也许我们会有一些重复的列）：

plot(s$d)

这意味着$V$的最后一列对于$Y$的重建非常小。为了说明这一点，我们可以考虑$V$最后一项为0的这种极端情况。在这个案例中，$V$的最后一列根本用不到。由于SVD的这种创建方式，$V$的列对$Y$的重建影响越来越小。我们通常认为这种描述为“解释了较少的变异”。这就意味着，对于一个大型矩阵，当你到达最后一列时，可能已经没有太多需要“解释”的内容了。例如，当我们把最后4列删除，看一下计算结果：

k <- ncol(U)-4
Yhat <- U[,1:k] %*% D[1:k,1:k] %*% t(V[,1:k])
resid <- Y - Yhat
max(abs(resid))

结果如下所示：

> max(abs(resid))
[1] 3.552714e-14

最大的残差基本上就等于0了，就意味着Yhat实际上是与Y一样，但是，我们至少需要4个维度来传输信息。

通过查看$d$，我们可以看到，在这个特定的数据集中，我们能得到一个很好的近似值，它只保留了94列。在下面的图形中，我们可以看到每列能解释的变异程度是多少：

plot(s$d^2/sum(s$d^2)*100,ylab="Percent variability explained")

还可以看一下累积曲线，如下所示：

plot(cumsum(s$d^2)/sum(s$d^2)*100,ylab="Percent variability explained",ylim=c(0,\
100),type="l")

虽然刚开始的时候，我们的数据是189维，但是我们可以使用95维来近似表示$Y$，如下所示：

k <- 95 ##out a possible 189
Yhat <- U[,1:k] %*% D[1:k,1:k] %*% t(V[,1:k])
resid <- Y - Yhat
boxplot(resid,ylim=quantile(Y,c(0.01,0.99)),range=0)

因此， 我们只使用了一半的维度就保留了原始数据中的大部分的变异：

var(as.vector(resid))/var(as.vector(Y))

结果如下所示：

> var(as.vector(resid))/var(as.vector(Y))
[1] 0.04076899

这个计算结果说明，我们使用了降维后的数据解释了原始95%的变异，我们需要注意的是，我们是通过$D$来计算的这个比例，如下所示：

1-sum(s$d[1:k]^2)/sum(s$d^2)

结果如下所示：

> 1-sum(s$d[1:k]^2)/sum(s$d^2)
[1] 0.04076899

因此，$D$中的元素可以告诉我们每个PC在解释变异方面所贡献的程度大小。

高度相关数据

为了辅助理解SVD是如何工作的，我们使用两组高度相关的列来构建一个数据集，如下所示：

m <- 100
n <- 2
x <- rnorm(m)
e <- rnorm(n*m,0,0.01)
Y <- cbind(x,x)+e
cor(Y)

结果如下所示：

> cor(Y)
          x         x
x 1.0000000 0.9998873
x 0.9998873 1.0000000

在这个案例中，第2列添加了很少的“信息”，因此所有的Y[,1]-Y[,2]都接近于0。使用rowMeans(Y)计算更加有效，这是因为Y[,1]-rowMeans(Y)和Y[,2]-rowMeans(Y)更接近于0。rowMenas(Y)最终生成的结果在$U$的第1列中。SVD的计算结果表明，仅使用第1列就能解释大多数的变异：

d <- svd(Y)$d
d[1]^2/sum(d^2)

结果如下所示：

> d <- svd(Y)$d
> d[1]^2/sum(d^2)
[1] 0.9999441

在这个案例中，许多列的数据高度相关，我们可以进行更大程度的降维操作：

m <- 100
n <- 25
x <- rnorm(m)
e <- rnorm(n*m,0,0.01)
Y <- replicate(n,x)+e
d <- svd(Y)$d
d[1]^2/sum(d^2)

结果如下所示：

> d[1]^2/sum(d^2)
[1] 0.9999065

练习

P338

投影

原始Rmarkdown文档参见作者的Github。

前面我们已经详细地描述了降维的概念，以及SVD和主成分分析的内容，现在我们来谈一下它们背后的数学原理。我们先从投影(projection)开始讲起。投影是一个线性代数的概念，它能帮助我们理解许多关于高通量数据的许多数学操作。如果想要了解更多相关的知识，可以找本线性代数的书来看一下有关投影的内容。在这一部分里，我们会提供一个快速的回顾，然后提供一些数据分析的相关案例。

作为回顾，我们需要注意的是，投影就是点与其子空间之间的距离，如下所示：

在上图中，顶部的点指向空间中的一点。在上图的这个卡通图中，空间是二维的，但是我们可以更加抽象地思考一下。这个空间由笛卡尔平面表示，小人站的这条线是点的一个子空间。将点投影到这个子空间上所对应的位置，就是这个子空间上这个位置距离原点最近的点。几何学告诉我们，我们可以通过从点到子空间一条垂线（虚线）来找到子空间上的这点。小人站在这个子空间上，这个人从原点走到投影点的位置时，就是这个点投影到子空间后的坐标。

为了扩展投影的概念，我们可以使用标准矩阵线性符号来说明这个点， $\vec{y} \in \mathbb{R}^N$是一个N维空间的点，$L \subset \mathbb{R}^N$是一个更小的子空间。

案例：当N=2

先看一个案例，$Y = \begin{pmatrix} 2 \\ 3\end{pmatrix}$，我们可以画出这个向量的图形：

library(rafalib)
mypar (1,1)
plot(c(0,4),c(0,4),xlab="Dimension 1",ylab="Dimension 2",type="n")
arrows(0,0,2,3,lwd=3)
text(2,3," Y",pos=4,cex=3)

我们可以马上定义一个坐标系统，将这个向量投影到空间中：$\begin{pmatrix} 1\\ 0\end{pmatrix}$ （x轴）和 $\begin{pmatrix} 0\\ 1\end{pmatrix}$ （y轴）。 $Y$ 向子空间的投影可以通过点2和3分别进行定义：

$$\begin{align*} Y &= \begin{pmatrix} 2 \\ 3\end{pmatrix} \\ &=2 \begin{pmatrix} 1\\ 0\end{pmatrix} + 3 \begin{pmatrix} 0\\ 1\end{pmatrix} \end{align*}$$

我们可以说 $2$ 和$3$是向量$Y$的坐标，$\begin{pmatrix} 1\\ 0\end{pmatrix} \mbox{and} \begin{pmatrix} 0\\1 \end{pmatrix}$ 是它的基。

现在我们定义一个新的子空间。红线（后面我们会画出这个图形）是一个子集(subset)$L$，它由满足 $c \vec{v}$ with $\vec{v}=\begin{pmatrix} 2& 1\end{pmatrix}^\top$的点构成。那么 $\vec{y}$ 在$L$上的投影就是$L$上最接近于 $\vec{y}$ 的点。因此我们需要找一个向量$c$，它是位于 $\vec{y}$ 和$c\vec{v}=(2c,c)$之间最小的距离。从线性代数的知识我们可知，这些点之间的距离正交于空间：

$$(\vec{y}-\hat{c}\vec{v}) \cdot \vec{v} = 0$$

上面公式也可以写为：

$$\vec{y}\cdot\vec{v} - \hat{c}\vec{v}\cdot\vec{v} = 0$$

即：

$$\hat{c} = \frac{\vec{y}\cdot\vec{v}} {\vec{v}\cdot\vec{v}}$$ $$\hat{c} = \frac{\vec{y}\cdot\vec{v}} {\vec{v}\cdot\vec{v}}$$

这里我们使用点号 $\cdot$ 来表示点积(dot product): $\,\, \vec{x} \cdot \vec{y} = x_1 y_1+\dots x_n y_n$.

下面我们使用R来演示一下上面的案例：

projection, fig.cap

mypar(1,1)
plot(c(0,4),c(0,4),xlab="Dimension 1",ylab="Dimension 2",type="n")
arrows(0,0,2,3,lwd=3)
abline(0,0.5,col="red",lwd=3) #if x=2c and y=c then slope is 0.5 (y=0.5x)
text(2,3," Y",pos=4,cex=3)
y=c(2,3)
x=c(2,1)
cc = crossprod(x,y)/crossprod(x)
segments(x[1]*cc,x[2]*cc,y[1],y[2],lty=2)
text(x[1]*cc,x[2]*cc,expression(hat(Y)),pos=4,cex=3)

我们需要注意的是，如果 $\vec{v}$ 满足 $\vec{v}\cdot \vec{v}=1$, 那么$\hat{c}$ 就是 $\vec{y} \cdot \vec{v}$ ，空间 $L$并没有发生改变，这处简化担任就是我们喜欢正交矩阵的一个原因。

案例：样本均值就是投影

设 $\vec{y} \in \mathbb{R}^N$ ，$L \subset \mathbb{R}^N$ 被以下向量张成：

$$\vec{v}=\begin{pmatrix} 1\\ \vdots \\ 1\end{pmatrix}; L = \{ c \vec{v}; c \in \mathbb{R}\}$$

在这个空间里，向量的分量(components）都是相同的数目，因此我们可以把这个空间看作为常数：在投影中，每个维度都是相同的值。那么$c$如何才能使得 $c\vec{v}$ 与 $\vec{y}$ 之间的距离最小呢？

当我们谈到这个问题时，我们会使用前面的二维图形。我们可以简单地抽象将$\vec{y}$视为N维空间上的一个点，将$L$视为一个更小数目的子空间，在这个案例中就是$c$。

回到我们的问题，我们知道，投影就是：

$$\hat{c} = \frac{\vec{y}\cdot\vec{v}} {\vec{v}\cdot\vec{v}}$$

在这个案例中它就是平均值：

$$\hat{c} = \frac{\vec{y}\cdot\vec{v}} {\vec{v}\cdot\vec{v}} = \frac{\sum_{i=1}^N Y_i}{\sum_{i=1}^N 1} = \bar{Y}$$

在这个案例中，它也非常容易使用微积分进行计算：

$$\frac{\partial}{\partial c}\sum_{i=1}^N (Y_i - c)^2 = 0 \implies 2 \sum_{i=1}^N (Y_i - \hat{c}) = 0 \implies$$ $$N c = \sum_{i=1}^N Y_i \implies \hat{c}=\bar{Y }$$

案例：回归也是一种投影

现在来看一下略微复杂的案例。简单线性回归也能用投影来解释。我们的数据 $\mathbf{Y}$（这里我们不再使用$\vec{y}$符号）是一个N维向量，我们使用一个线性方程$\beta_0 + \beta_1 X_i$来预测$Y_i$ 。此时我要找到能够使$Y$和由以下向量定义的空间的最小距离时的 $\beta_0$ 和 $\beta_1$ ，其中：

$$\vec{v}_0= \begin{pmatrix} 1\\ 1\\ \vdots \\ 1\\ \end{pmatrix} \mbox{ and } \vec{v}_1= \begin{pmatrix} X_{1}\\ X_{2}\\ \vdots \\ X_{N}\\ \end{pmatrix}$$

我们的 $N\times 2$ 矩阵 $\mathbf{X}$ 是 $[ \vec{v}_0 \,\, \vec{v}_1]$ ，$L$中的任何点都可以被写为 $X\vec{\beta}$.

正交投影的多维形式的方程为：

$$X^\top (\vec{y}-X\vec{\beta}) = 0$$

我们在之前看到过种形式：

$$X^\top X \hat{\beta}= X^\top \vec{y}$$ $$\hat{\beta}= (X^\top X)^{-1}X^\top \vec{y}$$

它向$L$的投影就是：

$$X (X^\top X)^{-1}X^\top \vec{y}$$

旋转

相关的Rmarkdown参见作者的Github。

与投影相关的一个最常见的应用就是坐标旋转(coordinate rotations)。在数据分析中，简单的旋转可以很好地对数据进行可视化和解释。我们将会介绍旋转背后的数学原理，并且给出一些简单的数据分析案例。

前面我们使用了下面的例子：

$$Y = \begin{pmatrix} 2 \\ 3 \end{pmatrix} = 2 \begin{pmatrix} 1\\ 0 \end{pmatrix} + 3 \begin{pmatrix} 0\\ 1 \end{pmatrix}$$

我们注意到，$Y$的坐标是$(2,3)$，现在我们使用如下的代码：

library(rafalib)
mypar()
plot(c(-2,4),c(-2,4),xlab="Dimension 1",ylab="Dimension 2",type="n",xaxt="n",yaxt="n",bty="n")
text(rep(0,6),c(c(-2,-1),c(1:4)),as.character(c(c(-2,-1),c(1:4))),pos=2)
text(c(c(-2,-1),c(1:4)),rep(0,6),as.character(c(c(-2,-1),c(1:4))),pos=1)
abline(v=0,h=0)
arrows(0,0,2,3,lwd=3)
segments(2,0,2,3,lty=2)
segments(0,3,2,3,lty=2)
text(2,3," Y",pos=4,cex=3)

但是，我们可以用其它的一些线性组合来表示点$(2,3)$：

$$\begin{align*} Y &= \begin{pmatrix} 2 \\ 3\end{pmatrix} \\ &= 2.5 \begin{pmatrix} 1\\ 1\end{pmatrix} + -1 \begin{pmatrix} \phantom{-}0.5\\ -0.5\end{pmatrix} \end{align*}$$

新的坐标就是：

$$Z = \begin{pmatrix} 2.5 \\ -1 \end{pmatrix}$$

从图形上我们可以看出来，这个坐标就是我们由新的基定义的空间的投影

library(rafalib)
mypar()
plot(c(-2,4),c(-2,4),xlab="Dimension 1",ylab="Dimension 2",type="n",xaxt="n",yaxt="n",bty="n")
text(rep(0,6),c(c(-2,-1),c(1:4)),as.character(c(c(-2,-1),c(1:4))),pos=2)
text(c(c(-2,-1),c(1:4)),rep(0,6),as.character(c(c(-2,-1),c(1:4))),pos=1)
abline(v=0,h=0)
abline(0,1,col="red")
abline(0,-1,col="red")
arrows(0,0,2,3,lwd=3)
y=c(2,3)
x1=c(1,1)##new basis
x2=c(0.5,-0.5)##new basis
c1 = crossprod(x1,y)/crossprod(x1)
c2 = crossprod(x2,y)/crossprod(x2)
segments(x1[1]*c1,x1[2]*c1,y[1],y[2],lty=2)
segments(x2[1]*c2,x2[2]*c2,y[1],y[2],lty=2)
text(2,3," Y",pos=4,cex=3)

我们可以使用矩阵乘法在表示$(2,3)$的这两个坐标中进行转换：

$$Y = AZ\\$$ $$A^{-1} Y = Z\\$$ $$A= \begin{pmatrix} 1& \phantom{-}0.5\\ 1 & -0.5\end{pmatrix} \implies A^{-1}= \begin{pmatrix} 0.5& 0.5 \\ 1 &-1\end{pmatrix}$$

其中，$Z$和$Y$表示了相同的信息，但是它们位于不同的坐标系中。

案例：双胞胎身高

我们先来看100个二维数据点$Y$，如下所示：

library(MASS)
n = 100
mypar()
set.seed(1)
y=t(mvrnorm(n,c(0,0),matrix(c(1,0.95,0.95,1),2,2)))
plot(y[1,],y[2,],xlab="Twin 1 (standardized height)",ylab="Twin 2 (standardized height)",xlim=c(-3,3),ylim=c(-3,3))

这里就使用了旋转：$Z = A^{-1} Y$，图形如下所示：

A = matrix(c(0.5,1,0.5,-1),2,2)
z = A%*%y
mypar()
plot(z[1,],z[2,],xlab="Average",ylab="Difference",xlim=c(-3,3),ylim=c(-3,3))

我们在这里进行的操作就是对数据进行旋转，从而使行新的坐标系$Z$的第一维是平均身高（就是相当于x轴），崦第二维则是两个双胞胎身高的差值（y轴）。

我们已经使用了奇异值分解计算主成分。有时候将SVD视为应动力非常有用，例如 $\mathbf{U}^\top \mathbf{Y}$ 就会构建出一个新的坐标系 $\mathbf{DV}^\top$ ，在这个新的坐标系中，它们的维度按照维度能够解释变异的程序进行排序。